classdef model    < handle
 properties
     samplingTime = 1;
    m_utilization, m_servers, m_num_users;
    price_per_vm_hour;
    vm_instantiation_time;
    
    %-------------
    d_num_users=[];
    mm=[];
    e1=[];
    e2=[];
    X;
    t=0;
    acc_pu=[];
    acc_inst=[];
    util=[];
    fired_inst = [];
    fired_t=[]; 
 end
 
 methods
     function obj=model
         obj.t=0;
         obj.samplingTime = 60;
         obj.price_per_vm_hour=0.08; %dollars
         obj.vm_instantiation_time=8;  %minutes         
     end
     
     function past_acc_inst=past_inst(obj)
         past_acc_inst = obj.acc_inst(:,1:size(obj.util,1))';      
     end
     
     function acc_inst=future_inst(obj)
         acc_inst = obj.acc_inst(:,size(obj.util,1)+1:end)';      
     end
     
     function past_acc_inst=past_future_inst(obj)
         past_acc_inst = obj.acc_inst';      
     end

     function past_pu=past_pu(obj)
        past_pu = obj.acc_pu(:,1:size(obj.util,1))';     
     end
     
     function past_util=past_util(obj)
          past_util = obj.util;     
     end
     
    function [m_x,d_x]=linearize(obj,x)
        m_x = mean(x);  
        d_x = x-m_x;  
    end

    function d_x=linearize_around(obj,x,m_x)
        d_x = x-m_x;  
    end

    % transfer function model of the behavior: workload+servers-> utilization
    % x=[num_users; num_servers]
    % y = avg_utilization
    % State-space model:  
    %         x(t+Ts) = A x(t) + B u(t) + K e(t)
    %         y(t) = C x(t) + D u(t) + e(t)
    function mm=identify_system(obj,data)        
        A=[1,-1; 0,1];
        B=[2,3;4,5];
        K=[1;1];
        C=[100,1];
        D=[100 1];

%           A=[1 -1 1 
%                0  1 2
%                0  1 2];
%          B=[2,3
%               4,5
%               4,5];
%          K=[4
%               5
%               5];
%          C=[1,1,1];
%          D=[1 1];
        
        m=idss(A,B,C,D,K,'Ts',obj.samplingTime,'DisturbanceModel','Estimate','InitialState','Estimate');
        %m1=init(m);
        %m.As=[NaN,NaN;NaN,NaN];
        % m.Bs=[NaN,3;NaN,5];
        % m.Ks=[NaN;5];
        % m.Cs=[NaN,NaN];
        % m.Ds=[NaN NaN];
        % m.x0s=[0; 0 ;1];
        m1=init(m);
        disp('building simulation model:.......');
        mm=pem(data,m1)  
    end
    
    function u=combineInputs(obj, d_servers_test , d_num_users)
        u =  [d_servers_test  d_num_users]; 
    end
    
    function mm=identifyModel(obj)
        load a_sample.mat
        intrvl=(1:219);
        sample_run_to_use=1;
        
        % linearization 
        utilization=utilization(sample_run_to_use,intrvl)' ;
        servers=servers(sample_run_to_use,intrvl)' ;
        num_users= num_users(intrvl)'; 
        [obj.m_utilization,d_utilization]= obj.linearize(utilization);
        %servers are also pu because  compute unit of 1 m1.smal instance (considered in data set) is 1
        % so we can assume we can pass compute units instead of  number of
        % servers
        [obj.m_servers, d_servers] = obj.linearize(servers); 
        [obj.m_num_users,d_num_users] = obj.linearize(num_users); 

        u =  [d_servers  d_num_users];
        y = d_utilization;
        % % u =  [servers  num_users];
        % % y = utilization;
        data=iddata(y,  u ,obj.samplingTime); 
        mm=obj.identify_system(data);        
    end
    
    function plot_compare_simple(obj,y1,y2)
        plot(obj.m_utilization + y1.OutputData,'g--');        hold on; 
        plot(obj.m_utilization + y2.OutputData,'r--');      
    end
    
    function plotit(obj,y,yn)
%         scatter(data.OutputData,y.OutputData);
%         plot( y.OutputData,'r--'); hold on; plot(data.OutputData);

        subplot(3,1,1); 
        plot(m_utilization + y.OutputData,'g--');        hold on; 
        plot(m_utilization + yn.OutputData,'r--');         hold on; 
        plot(m_utilization +data.OutputData);
        title('measured/modeled Utilizations')

        subplot(3,1,2); 
        plot(m_servers+data.InputData(:,1))
        title('#Servers')

        subplot(3,1,3); 
        plot(m_num_users+data.InputData(:,2))
        title('#Users')
    end

    % servers variable is a sequence over an interval of minute unit 
    % we are going to sweep through the time
    function servers=per_hour_cost_filter(obj,servers_) 
       % servers = [ones(1,45) 2.*ones(1,20)];
        servers=[0 servers_ zeros(1,60)];        
        T =size(servers , 2); 
         for t=1:T+60             
                if (t<T && servers(t+1)>servers(t)) %any time there is a scaleup
                    disp('hi')
                    % make sure the result of scaleup is accounted for, for 60 minutes
                    servers(t+1:t+60) = max(servers(t+1:t+60),servers(t+1)); 
                end
         end         
    end
    
    function cost_=cost(obj,servers_)
        hours=(sum(obj.per_hour_cost_filter(servers_))/60);
        cost_ = hours * obj.price_per_vm_hour;
    end
    
    % response time for a specific service rate
    function rt_=rt_static(obj,util,mu)
        rt_=(1/(1-util))*(1/mu);
    end    
    
    % this gives you almost the cdf
    % to get inverse sample for a population of users 
    % use 1-prob as a cdf 
    function prob_=prob_of_rt_more_than(obj,rt,T)
        prob_ = exp(1).^(-T/rt);
    end
    
    % this will take effect in 8 minutes
    function how_many_to_add=decide(obj)
        
    end
    
    function simuate(obj)
        intrvl=(1:219);
        obj.vm_instantiation_time;
    end
        
    % you decide on the number of servers 
    function util=getNextUtil(obj,server_num)
        % a hack to make pu s scale        
        % scaling down pu around linearization point
        % take [1,240] unit and map it to [1,6] unit
        % <0 is exception and >240 is not allowed
        % pu as constraint in algorithm
       ceiling_pu=240;
       tmp=server_num;
       server_num=1+(5*((server_num-1)/ ceiling_pu));
        
        
        m=obj.mm; 
        x=obj.X;        
        %A, B, C, D, K, and X0:
        %linearize this value around the mean of model because that the way
        % the model was built 
        d_server_num = (server_num - obj.m_servers); 
        % debug 
        % disp(strcat('d_server_num',num2str(d_server_num),'at',num2str(obj.t),' was ',num2str(tmp),' is ',num2str(server_num)));
        % disp(strcat('t:',num2str(obj.t)));
        u_t=[d_server_num ; obj.d_num_users(obj.t)];
        y_t = m.C * x + m.D * u_t + obj.e1(obj.t);        
        obj.X = m.A * x + m.B * u_t  + m.K*obj.e2(obj.t); %x(t+Ts) = A x(t) + B u(t) + K e(t)
        % de-linearize y_t
        util = obj.m_utilization + y_t;
        % hack for when model is not good for high capacity
%         if util<0
%             util=0;
%         end
         
        obj.util=[obj.util ; util];
        obj.t=obj.t + 1;              
    end
    
    function util=getNextUtilAcc(obj,instance_quantity_to_fire)
        a=instance_quantity_to_fire*prices_amz().get_compute_units;
        t=obj.t;
        if (t+60>size(obj.acc_pu,2))
           obj.acc_pu=[obj.acc_pu zeros(1,t+60-size(obj.acc_pu,2))];          
           obj.acc_inst=[obj.acc_inst zeros(size(instance_quantity_to_fire,2),t+60-size(obj.acc_inst,2))];          
        end
       obj.acc_pu(t:t+60)=obj.acc_pu(t:t+60)+a;
       obj.acc_inst(:,t:t+60)=obj.acc_inst(:,t:t+60)+instance_quantity_to_fire'*ones(1,61);
       if not(isempty(find(instance_quantity_to_fire>0)))
           obj.fired_inst = [obj.fired_inst; instance_quantity_to_fire]; 
           obj.fired_t=[obj.fired_t; obj.t]; 
       end
       % disp(strcat('instance quantity is ',num2str(obj.acc_pu(t))));
       util=obj.getNextUtil(obj.acc_pu(t));                  
    end
    
    function n=initialize(obj,length_in_minutes)
        obj.mm=obj.identifyModel;
        obj.t = 1;              
        obj.X=obj.mm.X0;
            
        % load the input of proper size and linearize it
       % length_in_minutes=219;         % assume 
      %  load a_sample.mat;
       % n= num_users((1:length_in_minutes))';          %load length_in_minutes of data from fifa
        n=workload().get_workload('data/day42_per_min.txt', 1, length_in_minutes);
        
        %generate random errors for the model
        randn('seed',0);
        obj.e1 = randn(length_in_minutes,1);
        obj.e2 = randn(length_in_minutes,1);
        obj.d_num_users = obj.linearize_around(n' , obj.m_num_users);        
    end

    function testGetNextUtilAcc(obj)
        duration=2*60;
        workload = obj.initialize(duration); 
         subplot(2,1,1);   plot(workload);
        utils=[];
        for i=1:duration-1
            if i==1 
                instance_quantity_to_fire=[0 4 0 0 0 0 0 0 0 0 0 ]; 
            else
                instance_quantity_to_fire=[0 0 0 0 0 0 0 0 0 0 0 ];             
            end
            % disp(strcat('time is ',num2str(i)));
            utils = [utils  obj.getNextUtilAcc(instance_quantity_to_fire)];      %one server
        end
        subplot(2,1,2);  plot(utils,'g--');          
    end
    
    function test_gen_sample(obj)  
        duration=23*60;
        workload = obj.initialize(duration); % the model has got identified here       , I want the simulator to initialize for 219 minutes of workload        
        subplot(3,1,1);   plot(workload);

        utils=[];
        for i=1:duration-1
           % disp(strcat('i:',num2str(i)));
            utils = [utils  obj.getNextUtil(4);];      %one server
        end
    
        % comparing with old one
%        intrvl=(1:duration);
        % mm=obj.identifyModel;
         load a_sample.mat;
         subplot(3,1,2);
        plot(num_users(1:219)');
        
        d_num_users = obj.linearize_around(workload , obj.m_num_users);        
        % modeled utilization with a new sequence of server num
        % d_servers_test2= obj.linearize_around(servers(1,intrvl)', obj.m_servers);
        d_servers_test2= obj.linearize_around(ones(duration,1)*4, obj.m_servers);
        u = obj.combineInputs(d_servers_test2,  d_num_users);
        % u
        data=iddata([],  u , obj.samplingTime);                         
        temp= sim(obj.mm,data,'noise');
        y_orig1 = obj.m_utilization+temp.OutputData; % simulates with input data
        subplot(3,1,3);
        plot(utils,'g--');  
        hold on; 
        plot(y_orig1,'r--');      
    end
    
    function mm=all(obj)
        intrvl=(1:219);
        mm=obj.identifyModel;
        
        % dont forget to indice it first sample_run_to_use=1;
        load a_sample.mat        
        d_num_users = obj.linearize_around(num_users(intrvl)' , obj.m_num_users);        

        % modeled utilization with a new sequence of server num
        d_servers_test1= obj.linearize_around(([3.*ones(1,100) 5.*ones(1,119)]   )',obj.m_servers);
        u = obj.combineInputs(d_servers_test1,  d_num_users);
        data=iddata([],  u , obj.samplingTime);                         
        y_orig1 = sim(mm,data); % simulates with input data
        
         % modeled utilization with original sequence of server num
        d_servers_test2= obj.linearize_around(servers(1,intrvl)', obj.m_servers);
        u = obj.combineInputs(d_servers_test2,  d_num_users);
        data=iddata([],  u , obj.samplingTime);                         
        % y_orig = sim(mm,data); % simulates with input data
         y_orig = sim(mm,data); % simulates with input data with noise       

        % original utilization from the dataset
        utilization=utilization(1,intrvl)' ;
        [obj.m_utilization,d_utilization]= obj.linearize(utilization);
        y = d_utilization;
        data=iddata(y,  [] ,obj.samplingTime); 

        obj.plot_compare_simple(data,y_orig);
        
        % yn = sim(mm,data,'noise'); % simulates with input data with noise       
    end    
    
    function plot_data(obj)
        load a_sample.mat        
        figure;
        subplot(3,1,1);   plot(num_users);
        subplot(3,1,2);  plot(servers(3,:)')
        subplot(3,1,3);  plot(utilization(3,:)')
    end
        
 end
end

